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Computational  Procedures  for  the 
Determination  of  a  Simple  Layer  Model 
of  the  Geopotential  from  Doppler  Observations 


Bertold  U.  Witte 

ABSTRACT.   The  geopotential  of  the  earth  is  rep- 
resented by  the  potential  of  a  simple  layer  dis- 
tributed over  the  earth's  surface.   Density  val- 
ues of  this  layer  for  104  surface  elements  have 
been  determined  from  Doppler  observations.   Com- 
putational procedures  using  this  method  are  pre- 
sented here.   Numerical  integration  is  employed 
to  compute  satellite  orbits  with  time  sequenced 
Doppler  observations  given  as  input.   The  geo- 
potential forces  acting  upon  the  orbit  are  de- 
scribed as  well  as  the  attraction  of  the  sun  and 
moon,  solar  radiation  pressure,  and  air  drag.   The 
results  of  the  orbit  fits  are  combined  in  an  ad- 
justment to  determine  the  density  values  of  the 
104  surface  elements. 


1 .  Introduction 

The  determination  of  the  earth's  gravity  field  from 
Doppler  observations  using  a  simple  layer  model  follows  an 
approach  outlined  by  Koch  [1968].   The  first  results  ob- 
tained by  this  method,  from  optical  satellite  observations, 
have  been  published  by  Koch  and  Morrison  [19  70].   These 
results  showed  good  agreement  with  existing  solutions ,  so 
this  approach  was  applied  to  other  available  data. 

Many  stations  of  the  worldwide  satellite  triangulation 
network  [Schmid,  1969]  are  close  to  the  Doppler  tracking 
sites  of  the  U.  S.  Navy  Doppler  Tracking  Network  (TRANET) 
[Anderle  1965]  and  connected  by  local  survey  to  the  BC-4 
stations  of  the  worldwide  triangulation  network.   In  order 
to  combine  the  results  of  both  systems  the  Doppler  observa- 
tions are  as  a  first  step  processed  by  means  of  the  simple 
layer  potential.   Results  are  given  by  Koch  and  Witte 
[19 70 ],  whereas  the  computational  procedures  are  described 
here . 

2 .  Observational  Equations 

Let  r  be  the  vector  of  the  apparent  position  of  the 
satellite  and 

r   =  [x,y,z]  (2-1) 

where  x,  y,  z  are  the  geocentric  rectangular  coordinates  of 


the  satellite;  the  origin  is  the  center  of  mass  of  the  earth; 
the  z-axis  is  identical  to  the  instantaneous  axis  of  the 
earth,  and  the  x-axis  points  at  an  angle  east  of  the  true 
vernal  equinox  which  equals  the  precession  and  nutation  in 
right  ascension  since  1950.   Since  the  orbit  computations 
extend  over  arcs  not  longer  than  7  days  this  coordinate  system 
is  a  good  approximation  to  an  inertial  system. 

The  rectangular  coordinates  of  the  station  vector 
(u  ,v  ,w  )  are  given  in  a  geocentric  coordinate  system  whose 

o     o     o 

orientation  coincides  with  the  worldwide  satellite  triangula- 
tion  network;  that  is,  the  w-axis  points  towards  the  mean 
pole  1900-1905  and  the  u-axis  towards  the  intersection  of 
the  Greenwich  meridian  (i.e.  the  zero  meridian  of  the 
Bureau  International  de  l'Heure  UTl-System)  with  the  equator. 

By  m<5ans  of  the  polar  motion  as  determined  by  the 
Bureau  de  l'Heure  and  the  sidereal  angle  as  defined  by 
the  Smithsonian  Institution  [1966]  the  station  coordinates 
are  rotated  into  the  x,y,z  system. 
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a 
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1 

(2-2) 


where  a,  b  are  the  components  of  the  polar  motion,  0  is  the 
sidereal  angle,  a  simple  linear  expression  of  UT1 ,  or 

0  =  0.rev277987616  +  1 .rev0027 3781191  (MJD-32 28 2 . 0 ) , 

MJD  is  the  Modified  Julian  Day  Number  as  defined  by  the 

Smithsonian  Institution. 

The  vector  r  of  the  apparent  position  at  a  time  t  is 

a  function  of  the  vector  e  of  the  orbital  elements  at  the 

epoch  t0,  the  vector  x  of  the  parameters  of  the  gravity 

field,  the  vector  of  the  station  coordinates  r  ,  the 

s 

radiation  pressure  parameter  KR  over  the  observation  interval 
and  the  drag  parameter  CL  over  the  same  time  interval. 

r    =    ?(e(t0),x,  ?s,  KR(t),  CD(t))  (2-3) 

For  the  formation  of  the  observation  equations  the 

following  expression  is  used,  which  is  derived  from  the 

Doppler  equation. 

f,   , 

f  =  fK  -  —  -Sr  ( |?  -  ?  | )  +  6f .  (2-4) 

b    c  dt   '      s  '       tro 

Here  f  is  the  calculated  frequency,  f,     the  base  frequency, 

c  the  velocity  of  light,  which  is  299,792,500  m/sec  [Fricke 

et  al.   1965],  -rr-lr  -  r  !  the  rate  of  change  of  the  distance 
'  dt '      s  '  & 

between  the  satellite  and  the  observation  station,  and  <$f 
the  tropospheric  refraction  correction.   A  detailed  de- 
scription of  the  model  used  for  tropospheric  refraction  is 
given  by  Witte  [1970],  but  a  condensed  version  is  included 
as  chapter  4. 


The  observed  Doppler  frequency  f0  is  obtained  by  comparing 
the  incoming  satellite  signal  with  a  reference  signal  generated 
from  a  local  precision  oscillator  and  measuring  the  time  required 
to  count  a  preset  number  of  beats.   In  order  to  get  the  calculated 
frequency  f(eq   (2-1))  the  base  frequency  f ,  ,  which  is  the  es- 
timated satellite  oscillator  frequency,  is  used.   The  assumed 
value  of  f,  may  contain  an  unknown  systematic  error  which  in 
this  report  is  considered  constant  over  one  pass.   The  observed 
frequencies  are  influenced  by  this  systematic  error,  which  is 
introduced  into  the  adjustment  as  a  bias  parameter. 

r.  ■*■     ■*■  "*■ 

Approximate  values  for  e,  x->    v    •>    K„ ,  C~  are  available  so 
that  we  obtain  by  means  of  a  Taylor  series  expansion  the  follow- 
ing observation  equation,  where  a  bias  parameter  b  for  the 
frequency  offset  of  every  base  frequency  for  each  pass  is  added 

Af.    =       I      (|£)        Aev      +    1E4   (||)         AX,    +       I      l^-)        Ar 

1         k=lV3e7ik         k°         £  =  lV^i£         l         m3lV8rs;. 


s 

mq 
lmq     ^ 


i     o        i 

i  denotes  here  the  ith  observation 

j  denotes  here  the  jth  pass 

q  denotes  here  the  qth  station  (total  of  n  stations) 

k  denotes  here  the  kth  orbital  element 

I  denotes  here  the  ilth  surface  density  element 

m  denotes  here  the  mth  coordinate  of  r 

g 

o  denotes  here  the  oth  orbit 

Af  =  observed  frequency  f0  minus  calculated  frequency  f. 

3 .   Determination  of  Partial  Derivatives  of  the  Doppler-Equation 
3.1   Evaluation  of  Derivatives  of  Doppler  Frequency  With 

Respect  to  Station  Coordinates. 

8  f 
For  the  partial  derivative   ~ using  equation  (2-4) 

sm 


9f       b   3(|r  -  rR  J  * )  (3-D 


Here 


9r        c        8r 
sm  sm 


d(  r  -  r   ) 
->    -*■    1  •  s  ' 

r  -  r 


s  '         dt 
and  the  slant  range  is 


I?  -  r  I  =  p  =  /(x  -  x  )2  +  (y  -  y  )2  +  (z  -  z  )2     (3-2) 

1      s  '  s       J         J  s  s 

The  differentiation  of  the  slant  range  with  respect  to  the 
time  yields 


d/Cx    -    x    )2    +    (y    -    y    )2    +    (z    -    z    )2 

s  J         J  s  s 


dt 
rtx  dv  dz  dx9  dyq  dzc= 

2/x-x    )2    +    (y-y    )2    +    (z-z    )2 
s  J    -'s  s 

(r   -   r    )    •     (r    -   r  , 

s    s ^_  (3-3) 
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1  s  ' 


Introducing  p  as  a  unit  vector  the  change  in  the  slant 
range  yields  with  equation  (3-3)  the  following  relation 


r  -  r  I"  =  p-(r  -  r  )  (3-4) 

s  '  s 


From  this  we  get  for  equation  (3-1) 
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with 
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where  co  is  the  angular  velocity  of  the  earth. 
Finally  we  get 
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3.2   Evaluation  of  Derivatives  of  Doppler  Frequency  with 
Respect  to  Parameters  Influencing  the  Orbit. 
The  derivatives  with  respect  to  the  parameters  e  ,  Xo»  ^ 

and  C  ,  as  used  in  equation  (2-5),  are  now  derived.   We  set 


3f 


3p 
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,  3X   3p 
=  1    m   r 


3X      3   a.   3X 
m     „   3f   __m 

m=l  3X   3pK 


(3-7) 


K 


m 


where  p„  stands  for  the  above  mentioned  parameters 

e^>  Xn?  Kt,  and  C~;  and  X   for  the  coordinates  x,  y,  z  of  the 
K    x,    R       D        m  '  J 

apparent  position  of  the  satellite.   The  partial  derivatives 

3X       3  X 

~ and  g  m  are  computed  with  the  help  of  variational  equations 

3pK      3pK  *  M 

(section  7  )  . 

The  derivation  of  the  Doppler  frequency  in  equation  (3-7) 
with  respect  to  the  apparent  position  of  the  satellite  leads  to 

r  D 
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Applying  (3-4) 
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Analogous  to  equation  (3-5)  and  with 
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3X 


0  we  get 
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The  differentiation  of  the  Doppler  frequency  with  respect  to 


the  velocity  X  ,  where  X   stands  for  the  coordinates  x,  y,z 
J      m  m 


leads  to 
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(3-14) 
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dr       9r 
For  the  components  of  the  derivative  -^—  and  ~ —  we  use 
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and 
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(3-16) 


3.3   Evaluation  of  the  Derivative  for  the  Bias  Parameter 
From  equation  (2-4)  we  get  for  the  partial  derivative  of  the 
Doppler  frequency  with  respect  to  the  base  frequency  f. 


3f  1  ,■>  -*■ 
ttp-  =  1  -  —  r  -  r 
df,         c  '       s 


(3-17) 


where  |r  -  r  |   is  obtained  from  equation  (3-3) 
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4 .   Tropospheric  Refraction  Model 

For  the  computation  of  the  correction  term  <5f     in 
equation  (2-4)  a  tropospheric  refraction  model  derived  by 
the  Naval  Weapons  Laboratory  is  used.   It  is  preferred  to 
models  developed  by  Hopfield  [1963  and  1969]  because  the  mean 
error  of  unit  weight  resulting  from  the  use  of  this  formula, 
which  will  be  derived  below,  is  somewhat  less  than  that  ob- 
tained with  the  Hopfield  models  [Witte  1971D 

The  range  error  As  caused  by  the  tropospheric  refraction 
may  be  described  by  the  following  equation 

As  =   Jn  ds  -  /  ds  (4-1) 

tro      g 

where  n  is  the  index  of  refraction  for  radio  waves  and  varies 
along  the  signal  path.   The  first  integral  is  taken  along  the 
signal  path  through  the  troposphere  and  the  second  along  the 
geometric  path.   If  both  ways  are  set  equal,  only  an  error 
of  second  order  occurs 

As  =  /  (n-1)  ds  (4-2) 

g 

The  contribution  of  the  tropospheric  refraction  to  the 
Doppler  shift  of  the  signal  received  from  a  passing  satellite 
is  given  by  (compare  also  Hopfield  1963) 

=  _  £b  d(As_).  (l+_3) 

tro      c   dt 
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In  equation  (4-3)  As  is  differentiated  with  respect  to  the 
time  t.  For  the  approximation  of  a  flat  earth  we  get  with 
equation  (4-2)  and  (4-3) 


H 
f       J(n-l)  dh 

5f    =  -  —  —  (° 

tro      c  dt      cos  Z 


) 


(4-4) 


The  angle  Z  is  the  zenith  distance  and  H  the  height  of  the 
satellite  above  the  earth.   After  integration  and  differentia- 


tion we  get  for  (4-4) 

f 


<Sf 


tro 


_b  (sin  Z   z  H  (-  _  1)) 
C       cos2Z 


(4-5) 


where  n  is  the  mean  refractive  index  and  Z  the  change  of  the 
zenith  distance  with  time.  In  order  to  solve  equation  (.4-5) 
the  following  relations  are  used  Csee  also  fig.  1) 

1)  The  component  of  velocity  along  the  direction  of 
the  unit  vector  q,  perpendicular  to  the  vector 

r  -  r   =  p,  and  in  the  plane  of  p  and  H,  is  q«p 

2)  The  component  of  velocity  along  q  can  also  be 

expressed  in  a  spherical  coordinate  system  centered 

at  the  observer  as   r  -  r  j  Z  or  pZ.   Thus 

i      s  i 


pZ  =  p.q 


Figure  1.   Diagram 

to  illustrate  eq-(4-6) 


(4-6) 


orbit 
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The  unit  vector  q  is  given  by  the  following  relation 

>x(pxH) 
px(pxH) 


q  =  PXrpXH)  (4-7) 


where 

|px(p  xH) |  =  sin  Z  (4-8) 

->- 

p  in  equation  (4-6)  can  be  rewritten  as 

->    -f        -f 

P  =  r  -  r  (4-9) 


and 


p  =  A  (u-io) 


Inserting  these  expressions  (4-7),  (4-8),  (4-9),  (4-10)  in 
equation  (4-6)  leads  to 


->- 


cos  Z         s     sin  Z 
and  substituting  px(pxH)  =    (p-H)p  -  (p-p)H  =  p  cos  Z  -  H 

into  (4-11),  equation  (4-5)  becomes 

f 

Sf.         =  -  -£  (r  -  r  )-(p L)  (n  -  1)     (4-12) 

tro      c        s    H    cos  Z 

We  replace  (n  -  1)  by 

.,     tro            AN,    -   N  .  , 
n   -  1  =  -4-  /      (n-l)dh  *  — ^- ^^       (4-13) 

H   h  .  .  (r  -  r  )-H 

1  '    stat  s 
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h  ,  .  =  height  of  station 
stat      & 

h^    =  height  of  troposphere 
tro      to         r 

AN     =  total  change  in  refractive  index  from  sea  level 

|! ';"'    up  to  the  top  of  the  troposphere  scaled  to  height. 

AN     =  change  in  refractive  index  from  sea  level  up  to 

the  height  of  the  station  scaled  to  height. 

AN  ,  .  is  determined  by  means  of  an  exponential  model 
s  ~ca~c 

similar  to  that  described  by  [Bean  et.  al .   1966]. 

C 

AN  .  .  =  C   (l-e"Clhstat)  +  — (4-14) 

stat     °  cos2Z 

C2 
In  order  to  get  flat  observations  corrected,  the  term  

cos2Z 

is  added.   Observations  with  Z>85°  are  deleted,  because 

the  computed  corrections  indicate  that  for  this  data  the  model 

used  is  overcorrecting. 

Combining  equation  (4-13)  and  (4-12)  leads  to  the  results 

given  in  the  following  expression 

f.  (   _  r-r     \  AN^_    -  AN  .  . 

of.    =  -  —  1  r  -  r    -  — [ = (4-15) 

tro      c   '      s  '    !■*■  i     r,   [     p  cos  Z 
I  r   cos  Z  J 

1  s  ' 

For  the  constants  C  ,  C, ,  C_  and  AN,    the  following  numerical 

o    1    2        tro  6 

values,  which  are  used  by  the  Naval  Weapons  Laboratory  (NWL), 

have  shown  good  results 

C   =  0.0025525 
o 

C1    =  0.1238 

C2  =  0.657X10"5 

AN^    =  0.0  02  3 
tro 

These  numbers  were  used  for  the  computation  of  the  tropospheric 

correction  model  as  given  by  equation  (4-15). 


"Private  communication  from  R.  J.  Anderle  [1970]. 
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5 .   Representation  of  the  Geopotential 

The  gravitational  potential  W  of  the  earth  is  divided 
into  the  normal  potential  U,  which  is  taken  as  known,  and 
the  disturbing  potential  T,  which  will  be  determined  by  the 
approach  used  by  Koch  [1968], 

W  =  U  +  T  (5-1) 

with 

7   n 


U  =  »i 

r 


1  +   E    E   C-)n  P    (sine}))  •  (  C    cos  mX 
n         n      r    nm  \  nm 

n=2  m=0 


+  S    sin  mX 
nm 


+  i  oo2  r2  cos2  <f>  (5-2) 


r,  <J> ,  A  are  spherical  coordinates  in  an  earth-fixed  coordinate 

system,  k  is  the  gravitational  constant,  M  the  mass  of  the 

earth,  a  the  mean  equatorial  radius,  P    the  fully  normalized 
'  H  '   nm         J 

associated  Legendre  function  of  degree  m  and  order  n,  and 

co  the  angular  velocity  of  the  earth.   C    and  S    are  fully 

nm      nm 

normalized  harmonic  coefficients.   Their  values  are  taken 
from  the  solution  of  Anderle  [1965].   The  coefficients  C2  1 
and  S2,i  have  been  deleted  because  of  the  chosen  coordinate 
system,  defined  above.   The  centrifugal  term  in  (5-2)  is  omitted 
for  points  outside  the  earth. 

The  disturbing  potential  in  equation  (5-1)  is  represented 
as  the  potential  of  a  simple  layer  with  the  density  x  C  <|) ,  X ) 
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distributed  over  the  surface  E  of  the  earth 

T  -_    jj   Xiiiil  dE  (5-3) 

E     S 

s  being  the  distance  between  the  point  at  which  T  is  com- 
puted and  a  variable  point  on  E.   To  evaluate  this  integral, 
the  earth's  surface  is  divided  into  elements  on  which  the 
density  of  the  simple  layer  is  assumed  to  be  constant  so  that 
the  integral  in  (5-3)  can  be  replaced  by  a  summation.   104 
surface  elements  AE,  are  chosen  here.   These  are  bordered  by 
parallels  and  meridians  and  approximate  the  size  of  a  20°  x  20 
area  at  the  equator. 

104  AT 

T  =  I      Xo   //   ~  (5-4) 

£=1      AE£  I 

The  integral  over  the  surface  elements  AEp  in  (5-4)  is 
computed  numerically  by  splitting  AEp  into  4  subdivisions 
and  assuming  the  kernel  of  the  integral  computed  for  the 
midpoint  of  the  subdivision  to  be  constant  over  the  subdivision 
To  compute  the  distance  s  between  a  satellite  and  the  midpoint 
of  the  subdivision  of  AE.  it  is  necessary  to  determine  the 
coordinates  of  these  midpoints.   For  this  purpose  the  equi- 
potential  surface  U(eq  5-2)  is  defined  as  an  ellipsoid  U0 . 
U0  is  computed  from  the  equation  which  holds  for  the  surface 
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of  a  level  ellipsoid 

kM 


U, 


f    2    ,2,1/2 
(a   -  b  ) 


arc  tan 


(a2  -  b2)1/2 


+  |  co2  a2   (5-5) 


All  the  parameters  of  (5-5)  are  available  from  (5-2) 
except  b.   We  obtain  b  from  b  =  a(l-f),  and  f  from  the 
value  of  C20  in  (5-2)  in  the  following  manner.   The  relation 
between  C2o  and  f  is  given  by  [Heiskanen  and  Moritz   1967, 
p.  110]. 

C20  =  -  J2 / *$ 

C5-6) 


t    1   e*2 
J2  ~-    3 


1+e 


'2 


,    _2_  me  ' 


15  qD 


with 


and 


1 
qQ  =  2 


1  +  — rsrj  arc  "tan  e 


■1 


m  = 


2  2 
(0  a 


kM/l+e'2 

'2  -    1 
"  (1-f)2 


-  1 


e'  is  obtained  from  (5-6)  as  a  function  of  J2  by  the 


Newton-Raphson  method  of  successive  approximation.   For 

dJ: 
de 


this,  the  derivative  -3 — f  is  needed 


dJ2  1    /  >       „'K     '),      ■/< 

de'    "    3  l 


2e'  3 
nrv2 


2 


T    2 


15      qoAi  +  e.2  Cl  +  e'34)^"    45    IT?7"2^  q0    kMv/(1  +  e,2)3 


m 


e'2         a)2a3 


me 


2q, 


((-i^)r^ 


6  4.  1    x       3 

pr-3-  arc   tan    e?    +   p-jr 


(5-7) 
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Now  all  parameters  for  the  computation  of  U0  ( eq   (5-5))  are 
known  and  therefore  the  radius  vector  of  the  reference  sur- 
face can  be  computed  in  setting  U0  =  U.   Solving  equation  (5-2) 
for  r  we  get 


r  = 


kM 


7 


n 


1  +   Z   Z      (-)n  P    (sin*)  x 

_0_nr    nm      T  ■*■ 
n=2  m=0 


(C    cos  mX   +  S    sin  mX) 
nm  nm 


co2r2  cos2(j) 
2U. 


(5-8) 


Equation  (5-8)  is  determined  for  the  earth-fixed  co- 
ordinates of  each  midpoint  by  direct  iteration  for  the  reference 
surface,  where  the  initial  value  for  r  is  obtained  for  a  given 
latitude  (J)  on  an  ellipsoid  of  revolution  with  semi-axes  a  and  b. 
To  these  values  of  the  radius  vector  the  topographic  heights 
above  sea  level,  published  by  Kaula  et  al .  [1966],  were  added. 
Employing  the  given  values  <f> ,  X   of  the  midpoints,  we  can  compute 
the  corresponding  u,  v,  w  coordinates,  and  then  the  x,  y,  z 
coordinates  of  the  midpoints  by  means  of  the  rotation  matrix  (2-2) 

6 .   Partial  Derivatives  of  the  Potential  W 

The  derivatives  of  the  force  acting  upon  the  satellite 
with  respect  to  the  position  of  the  satellite  are  needed  in 
the  solution  of  the  variational  equations  (Sec.  7)  for  the 
orbit  computation.   These  derivatives  are 


82W    32W     82W     82W    82W 


3x2  '  3x3y  '  3x3z  '  By7  '  8y3z 


,  and 


32W 


The  gradient  of  the  potential  U  is  used  for  the  computation 
of  the  reference  orbit.  The  approximate  density  values  are 
set  equal  to  zero.  Therefore  the  gradient  of  the  potential 
T  will  not  be  computed  here. 

6.1   Gradient  of  the  Potential  U 

Using  the  chain  rule  method  for  the  partial  derivatives 
of  U,  where  U  =  U(r,  cf>jA),  we  find 


3x    3r  3x    8cj)  3x    3X  3x 


(6-1) 


M     £uto+3U8£3U8x 

9y  '  8r  8y    8(f)  3y    9X  8y 


(6-2) 


m    min  +  mii  +  MiA 

3z    '     3r    3z         3<J>    3z         3X    3z 


(6-3) 


The  differentiation  of  equation  (5-2)  with  respect  to  r  yields 


3U 
3r 


kM  r 

"  "  r^L 


7  n 

1  +   Z  Z   (n+1)  (-)   P    (sin<J>)  x 

o n  r     nm 

n=2  m=0 


(C    cos  mX  +  S    sin  mX ) 

nm  nm 


+  w2  r  cos2cJ>  (6-4) 


d  P    (sincf>) 


With 


nm 


d(J> 


=  -m  tan  d>  P    (sin  d>)  +  P    ,,  (sin  <j>) 
Y   nm      Y     n,m+l 
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we  obtain 


1  r\ 

|x  =  -  —  [  E    2   C-)   (m  tan  cj>  P    (sin  cf>)  -  P    ...  (sin  <(>)) 
^      rLn=2  m=0   r  nm  n>m+1 

(C    cos  mX  +  S    sin  mX )  -  w2  r2  sin  d>  cos  d>      (6-5) 

nm  nm        J 

and 

£E  =  -  — [  E    S   m  (-)n  P    (sin  cf>)(C    sin  mA  -  S    cos  mX) 
8X     r    o    „    r    nm      T   nm  nm        J 

Ln=2  m=0 

(6-6) 
As  mentioned  earlier  the  centrifugal  terms  in  C6-4)  and 
(6-5)'  should  be  omitted  for  points  outside  the  earth. 
The  other  partial  derivatives — equations  C6-1)  -  C6-3)--are 


found  with  r  =  A2+y2  +  z2,  tan  <£  =   .      =  ,  and 

/x2  +  y2 


tan  X  =  J 


8x  " '  r*   8y  "  rs   3z  "  r  ; 


M         xz       8jb_  =     -yz      8_£  .  /x2  +  y2      ,.  0. 
3*  "  "  Wx2  J  y2  '  K    "  Wx2  +  y2  '  8z  '     r*         C    ^ 


|4  =  — ^ ,  |i  =  ^ ,  |i  =  0  (6-9) 

3x   x2  +  y2   9y   x2  +  y2   8z 


6.2   Second  Derivatives  of  the  Potential  U 

The  second  derivatives  with  respect  to  x,  y  z  of  the 
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potential  U  yield   (using  equations  (6-1)  -  (6-3))  the 
following  expressions 

32U    32U  /9rV 


9x- 


3r2 
+  2 


3x 


3ZU  /8(f) 


3<}> 


3x 


8_nj 

3A2 


3x 


32U   3r  3_£  +     32U   3r  3A_  +     32U   _3J>_  3_A 


3r •  3<j)  3x  3x 


3r«  3 A  3x  3x 


3<J>  •  3  A  3x  3x 


+  1H  ilz  +  1M  ill  +  8U  ilA 


3r    -    2 
3x 


3  9         rs 

3x 


3  A    _ 
3x 


C6-10) 


32U         3.2U    3r    3r         3^U    3^   3^         3JU    3jV_  _3J. 
3x3y         j,    2    3x    3y         3c})2    3x    3y    '     -.2    3x    3y 


32U 
3r3(J> 


32U 
3A3c}) 


/3r   ii  +    3r    3£\  +    3  2U    /3r   _3_A    +    3r    3_A\ 
\3x    3y        3y    3x/      3r3A\Jx    3y    '     3y    3x  / 

/  3A_  ii  +   il   d±\ 
\  3x    3y        3y    3x  / 


+    3U    32r      +    3U    32(j)      +    3U    32A 
3r    3x3y         3$    3x3y         3A    3x3y 


(6-11) 


32U    .     3^U    3r    3r        _3_^U    3£   3£   +    _3_fU    3_A    3_A 


3x3z 


3r 


2    3x    3z 


3<1) 


2    3x    3z 


3A 


2    3x    3z 


32U 
3r3(}) 


32U 
3A3(j) 


/  3r  2i  +  9r   3  $  \        32U     /3r    W    ,     3r    3_A  \ 

^ 3x    3z  3z    3x  /        3r3A^3x    3z         3z    3x  ) 

(  3_A    d±    .  _3A    Z±\ 

\3x    3z  3z    3x/ 


,     3U    32r       .     3U    324>       ,     3U    3  2A       f6-12) 
3r  "SxTz         3^  ?x3~z        "3~X  "3~xTz 


32U 


3y: 


32U  /  3rX2 


3r- 


3y 


32U 
3<})2 


ii 

3y 


32U 


3A 


3A^2 
3y, 


+  0  92u    lz  ii  +       92u     9r  li  +  0  32u    2*  ii 

3r3(J)    3y    3y  3r3A    3y    3y  3X3$    3y   3y 
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+  £U  lf£  +  £U  3J^_  +  3U  3J\X 

"  "  '  H    3y2 


3r  3y2 


3X  3y2 


C6-13) 


32U   _  ^U  3r  3r    SJ^U  3f  3|  +  S^U  3X.  3X 


3y3z  "  3r2  3y  3z      2  3y  3z 


3X 


2    3y    3z 


^!£_  f i£  i£   +    9r    3^\         32U     /  3r    3A     ,     3r    3X\ 
3r3(J>   \3y    3z         9z    3y  /         3r3X  \  3y    3z         3z    3yJ 


32U 


( IX    d±    .     d\    d±\    i  3U  32r   +  3U  32({)   +  3U  32X    (6-14) 
9A9<j>  \   3y  3z    3z  3yJ    3r  3y3z    3$  3y3z    3X  3y3z 


3ZU    32U 


3z; 


3r' 


+  2 


32U   3r  3 


3r3<£  3z  3z 


^  +  2 


32U   3r  3X 


~  ^  +  2 


3r3X  3z  3z 


32U   3_X  3^ 
3X3(b  3z  3z 


3U  32r    3U  3^1    3JJ  32X 

3r  3z2    3d)  .  2  '  3X  _  2 
T  3z        3z 


(6-15) 


Now  the  partial  derivatives  in  equations  (6-10)  -  (6-15) 
are  derived.   With  equation  (6-4)  we  find 


32U  ,  2kM 
3r2     r3 


1 


n 


1  -:  ;,       Z  (n  +  l)(n+2)  (-)"  P   (sin  4) 
2    0    -  r     nm      r 

n=2  m=0 


(C    cos  mX  +  S    sin  mX ) 
nm  nm 


(6-16) 


I^U   =  _  kM         (n+1)  (a5  (p-      (sin  4,) 

3r3d;       2    o   -n       r    n,m+l 
T      r  Ln=  2  m=0 


m  tan  d>  P   (sin  6)  ) 

T  nm      r 


(C    cos  mX  +  S    sin  mX ) 
nm  nm 


(6-17) 
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9r9X 


92U       vm  r  7    n  ^n- 

-  SI    E    E   (n  +  1)  (|)   P^   (sin<f))(-C    sin  mA 
r2  ln=2  m=0        r    nm  nm 


+  S    cos  mX) 
nm 


] 


(6-18) 


8ZU  __  kM 

,,2  "   r 


r  7    n 


n 


_n=2  m=0      \  cos  cf> 


E    E   (-)"(  -  -  --  P    (sin  cj))  +  m  tan  <j> 


(P      m+1    (sin    (J))    -   m   tan    <f>    P         (sin    <j>))      +    (m+1)    tan    <J> 

P«    mil     (sin    <f>)    "   ^      xo    ^sin    <J>))CC        cos   mX 
n,m+l  T  n,m+2  r J     nm 


+    S        sin   mX) 
nm 


(6-19) 


92U 
9<J)9X 


kM 
r 


r    7  n  n 

E  Em    (  — )       (m   tan    d>    P         (sin    <J>) 
o         n  2?  nm  v 

_n=2  m=0 


-  P    ,,  (sin  d)))(-C    sin  mX  +  S    cos  mX) 
n,m+l       r      nm  nm 


(6-20) 


32U 
9X2 


kM 
r 


"  7    n        n 

E    E   m2  (§)   P_  (sin  <|>)(C,._ 
-n=2  m=0 


r    nm 


nm 


+  S    sin  mX ) 
nm 


] 


(6-21) 


And  for  the  other  terms  in  equations  (6-10)  -  (6-15)  we  get 


9  2r  _  1    x2  .   9  2r  _    xy  .   9  2r  _    xz 

„  2    r     3  '  9x9y       3  '  9x9z       3 
9x        v  J  r  r 


(6-22) 


9  0 


9x2    r2/x2+y 


=  /_  1  +  2xi+_xf_\ 
+v2  V       r2    x2+v2  ' 


r>z   xz+y 
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92(j)   _    xyz    /  ,2  +    1   \     82(J)  _     -x  +    2xz2     (6-23) 

3x3y     r2/x2+y2  \T2  x2+y2)      '   9X8Z     r2/x2+y2     T» /^T^T 

^  -       2^        ;     *ii  =      y2-*2     ;     ^fA  =  0  (6_24) 

9x2         (x2+y2)2         dxdy         (x2+y2)2         9x8z 


92r        1        y^        92r  yz 

^2        ^  3    '    9y9z  3 

9y  r  J  r 


924>    :  z  /  2y2     _    .     +         y2 


(6-25) 


2         r2/? 


=  /2Yl  _  1  +  _r_\ 

+  v2  V     r2  x2+v2  / 


8y^        r"/x"+y"  V    r2  xz+y 


92<J>    .    _  Y  +         2yz2 

8y3z  r2/x2+y2        r*^1^ 


(6-26) 


ill  ■ 2xy__  .  _^A  =  Q  (6_27) 


^2      /  2 i  2\2    9y  9  z 
9y       (x  +y  ) 


92r    1  _  zj_    9j^_   _  2/x2+y2  z     92A  _  (6-28) 

9z2    r   r3  '  9z2         r"        8z2 


This  method  of  computing  the  partials  of  U  was  found  by  Gulick 
[1970]  to  be  the  most  efficient  compared  to  two  other  alternatives. 
Subroutines  which  compute  these  partials  based  on  this  and 
other  methods  may  be  found  in  this  reference  along  with  a  de- 
tailed discussion  of  comparative  efficiencies  and  a  complete 
derivation  of  expressions  given  in  equations  (6-10)  -  (6-28). 
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7 ,   Variational  Equations 

The  position  of  a  satellite  can  be  described  as  a 

function  of  the  parameters  of  the  earth's  gravity  field  plus 

non-conservative  forces  F 

c 

r  =  VW(r,x,t)  +  F  (r,r)  (7-1) 

For  the  time  being  the  non-conservative  forces  F   will  not 

&  c 

be  considered  here  so  that  we  have 

r  =  VW(r,x,t)  (7-2) 

For  the  solution  of  this  differential  equation  a  numerical 

integration  procedure,  as  described  in  the  next  section  is 

needed.   In  order  to  do  an  orbit  adjustment,  additional 

equations,  called  variational  equations,  are  necessary  for 

3X         9X 

the  determination  of  the  partial  derivatives   ~ and   ~ 

9pk        8pK 

given  in  equation  (3-7).   To  make  the  notation  more  concise, 


*8t  =  (?,?)  and  r0st  =  Cr0J?0, 


we  form  state  vectors  r__,_  =  (r,r)  and  r0    =  (r0  ,r0 )  .   r 


consists  of  the  quantities  X  ,  X  .   rn    constitutes  a  subset 

n  mm     °st 

of  the  parameters  p,  .   Here  r0  is  the  initial  position  vector 

?  ...  ->      -f 

and  r0  the  vector  of  the  initial  velocity;  r0  and  r0  are 

equivalent  to  the  six  orbital  elements  represented  by  e,  in 
equation  (2-5).   Using  the  state  vector  notation  we  get  for  eq  (7-2) 

?st  =  (r\VW)  =  S8(r8t,x»t)  (7-3) 


25 


The  position  of  a  satellite  is  now  given  as  a  solution  to 

the  differential  equation  (7-3)  with  the  initial  state 

vector  as  boundary  condition.   If  we  form  the  time  derivative 

of  the  partial  of  r  .  with  respect  to  the  parameters  p. 
r  St  -^  K 

we  get  J   , 

d  [~5l 

V%  1  ^st  (7-4) 

dt   "  .-*• 

The  numerical  integration  procedure  will  be  applied  to  (7-4) 
in  order  to  obtain  the  required  partial  derivatives.   First 
we  will  determine  the  above  mentioned  time  derivative  of 
the  partial  of  r   with  respect  to  the  initial  state  vector 
and  get,  with  eq   (7-3), 

V3r°2i/=  8p»t  3rst  ,  3*st  3x   ,  3*st   at 


dt        a  +     -HJ         +   --*  9t 

3r  .  9r„  ,     3y   9r0  .         dr0  , 
st    °st     A     °st  °st 


Note  that  the  last  two  terms  are  zero,  because  x  an^  t 
do  not  depend  on  the  initial  state  vector.   So  we  finally 
get 

/9?st  \ 

\  ost/_  =    st  st  (y_5) 


dt 

St      °St 


3r_   9r0, 
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Using  the  rectangular  coordinate  system  x,  y,  z,  defined  in 
section  2,  we  get  in  matrix  notation  for  the  first  factor 
in  eq   (7-5) 


9r 


st 


9r 


st 


9x 

9x 

9x 

9x 

9x 

9x 

9x 

9y 

9z 

9x 

9y 

9z 

9? 
9x 

9y 

9y 

9y 
9z 

9y 
9x 

9y 
9y 

9y 

9z 

9z 

9z 

9z 

9z 

9z 

9z 

9x 

9y 

9z 

9x 

9y 

9z 

9x 

9x 

9x 

95^ 

9x 

9x 

9x 

9y 

9z 

9x 

9y 

9z 

3£ 

9x 

9y" 
9y 

9z 

12 

9x 

9£ 
9y 

9£ 
9z 

9z 

9z 

9z 

9z 

9z 

9z 

9x 

9y 

9z 

9x 

9y 

9z 

(7-6) 


If  we  partition  this  matrix,  we  will  get  the  following 
structure 


9r 


st 


9r 


st 


0 


92W 


9X 


m 


9X.  9X.   ; 


^ 


(7-7) 


where  0  is  a  3^3  zero  matrix  and  I  is  a  3x3  identity  matrix 

-*■     "*■ 
(since  r  and  r  are  the  independent  parameters  of  the  differential 

equation  (7- 3)), the  lower  left  hand  portion  arises  from 

differentiation  of  (7-5)   with   92W   as  a  shorthand  notation 

9X.9X. 
i   ] 
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for  this  part,  and  the  lower  right  hand  portion  is  0  if  drag 

is  not  present  or  is  neglected;  terms  from  drag  would  occur 

in  the  entire  lower  half  of  the  matrix  if  drag  were  present. 

If  a  simple  radiation  pressure  model  is  included,  a  Dirac  delta 

function  term  occurs  in  the  lower  left  portion  of  the  matrix; 

for  a  very  complex  model  with  re-radiation  an  impulselike 

component  would  still  appear. 

In  setting  up  the  variational  equations  for  the  solution 

of  x  an<3  r0   ,  the  nonconservative  forces  F  are  not  being 
st  -"         <^ 


9X 


applied.   Hence,  we  set  the  term  ^  £n  equation  (7-7) 


9X 


m 


to  zero,  and  get  instead  of  (7-7) 


- 

3rst 

0 

I 

3?  . 
st 

92W 

0 

9X.9X. 

(7-8) 


The  differentiation  of  the  state  vector  r  .  with  respect  to 

st         r 

the  initial  state  vector  r0    yields 

st 


9x 


9x, 


dx 


ax 


9x 


9y< 


9z, 


dx, 


dx 


9yc 


9x 


9x0 

9yQ 

9z0 

9x0 

3yc 

9z0 

9y_ 
9x0 

9y 

ay0 

9y 
9z0 

9y 
9x0 

9y 

3yQ 

9y 

9z0 

3?  t 

st 

9z 
dx0 

9z 

9y0 

9z 
9z0 

9z 

ax0 

9z 

3yc 

9z 
9z0 

st 

dx 
9x0 

dx 

3yQ 

9x 
9z0 

9x 
9x0 

9x 

8yG 

9x 
9z0 

9y 
9x0 

3y 
9yD 

9y 
9z0 

3y 

3x0 

9y 

3y0 

3y 

3z0 

9z 

d± 

9z 

9z 

9z 

9z 

9z, 


(7-9) 
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For  the  start  of  the  integration  procedure  this  matrix  has 
the  following  structure 


3r 


st 


3?, 


st 


0 


(7-10) 


For  the  time  change  of  the  derivatives  of  the  state  vector  r 
with  respect  to  x  we  get,  analogous  to  (7-5), 


st 


3r 


st 


X 


dt 


^st  9?st 


3? 


3r 


st 


st 


ax 


ax 


(7-11) 


with 


3? 


st 


3X 


3x 


3X 


3x 


9-Xi 

3X2 

3y 

3Xi 

3y 
3X2 

3z 

3z 

3Xi 

3X2 

3x 

3x 

3Xi 

3X2 

By 
3Xi 

3y 

3X2 

3z 

dz 

3X2 


3x 


3Xi  04 

3Xi  04 

3z 
3Xi  o  k 

3x 


3X 


10  4 


*± 


3Xi  04 

3z 
3Xi  04 


(7-12) 
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3r 


st 


w 


o 

0 

0 

3x 
3Xi 

3y 

3Xi 
3z 

3Xi 


0    .  . 

.  .    0 

0    .  . 

.  .    0 

0    .  . 

.  .    0 

3x 

3X2   '  * 

3x 
3Xi  ot 

ay 

3X2   '  • 

3'y 

3Xi  o  4 

3z 

9z 

3X 


3X 


1  0  h 


(7-13) 


3r 


and 


st 


supplied  by  (7-8 ) 


3r 


st 


3r 


In  (7-11), 


st 


3X 


represents  the  explicit  instantaneous  de- 


pendence of  velocity  and  acceleration  on  the  gravity  field.   Since 

instantaneous  velocity  is  an  independent  variable,  representing, 

with  the  instantaneous,  position   the  instantaneous  state  of  the 

satellite,  the  top  half  of  (7-13)  is  null.   On  the  other  hand, 
3r 


st 


3X 


represents  the  dependence  of  the  state  vector  (position  and 


velocity)  on  the  gravity  field,  considering  the  state  vector  as 

obtained  by  integrating  the  accleration  over  time;  hence  the  terms 

of  (7-12)  are  not  necessarily  zero,  although  the  lower  half  of 

(7-12)  literally  appears  to  be  the  same  as  the  upper  half  of  (7-13). 

Initial  conditions  for  the  integration  are    st i   =  I  0      (7-14) 

o    L  -i 

For  the  elements  of  matrix  (7-13)  we  find 


3X, 


3x 
3X£ 

3y 

3X 

3z 

3X0 


c3?) 


3X0   3* 


•I 


3 


cli) 


3x0  3y 


■i 


(P) 


3X^   3z 


(7-15) 
(7-16) 
(7-17) 
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where  T  is  the  disturbing  potential  as  defined  in  eq.  (5-3). 


Using  the  distance  s.  =  Ax  -  x^)2  +  (y  -  yJ2  +  (z  -  z.)2 

between  the  satellite  and  the  midpoint  of  the  £th  surface 
element  we  get  for  the  partial  derivatives  of  T  with  respect 
to  x,  y,  z  from  eq   (5-4) 


104 

S 
£=1      AE, 


|£  =  I      X,      !!      "  ~^dE  (7-18) 


9T 


'I 
104 


*~xz 

dE 

s3 

y-y* 

dE 

s3 

Z"z£ 

dE 

s3 

3v  S   X£  //   -  -5T*  dE  C7-19) 

dy  £  =  1   *  AE£ 

~T  104 

f^  S   Xo  //   "  "TT1  dE  (7-20) 

Z  £=1   *  AE£ 


For  the  equations  (7-15)  -  (7-17)  we  now  get 


x-x. 


af  (^}  =  //  "  ~^dE  (7"21) 

dX£   dX      AE£     s 

if-  (H)  =    //    "  ^  dE  (7"22) 

dX£   dy      AE£     S 

-JL  (II)  =   // jA  dE  (7-23) 

d*£   dZ      AE£ 
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8 .   Orbit  Integration  Procedures 

The  equations  of  motion  and  the  variational  equations 
of  the  satellite  expressed  in  the  rectangular  coordinate  system 
defined  in  sec.  2  are  integrated  numerically  with  an  0.8  minute 
time  step.   The  twelfth-order  Cowell-Stormer  multistep  process 
is  used  for  the  positions  and  the  tenth-order  predictor 
corrector  formulas  of  Adams-Bashforth  and  Adams-Moulton  re- 
spectively are  used  for  the  velocities  [Henrici   1962,  pp.  191- 
198,  291-294],   The  accompanying  flow-chart  shows  the  struc- 
ture of  the  above  mentioned  procedures. 

The  derivatives  with  respect  to  the  initial  state 
vector  and  to  the  unknown  density  values  (i.e.?  eq   (7-5)  and 
(7-11))  are  obtained  numerically  by  integrating  the  variational 
equations.   In  this  integration  the  corrector  is  solved  ex- 
plicitly without  predicting,  using  the  linearity  of  the 
variational  equations  [Riley   et  al.   1967]. 

These  procedures  require  the  use  of  a  starting  process, 
because  the  first  time  steps  cannot  be  computed  without  knowing 
previous  values.   Gill's  modification  of  the  Runge-Kutta  in- 
tegration method  [Romanelli   196  0]  is  therefore  employed  for 
the  computation  at  these  time  steps.   The  methods  mentioned 
are  now  explained. 
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PREDICT  ?p+i 
ADAMS- BASHPORTH 


PREDICT  ?p+1 
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1    =   0 
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ADAMS -MOULTON 

?p+l  ?p+l 
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COMPUTE 
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Vp+1  SETS 
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COWELL 

yPp+1     Tp+1 

A>+1     Vp  +  1, 


Figure  2.  Integration  process  for  the  computation  of 
position,  velocity,  and  acceleration 
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8.1   Runge-Kutta-Gill  Method 

The  totality  of  the  variational  equations  is,  as  has 
been  shown  in  section  7,  a  6x(6+104)  matrix  with  a=aCr0  . ,x), 
where  a  is  a  subset  of  parameters  of  p. 


3(rst> 


8(?st) 


8(a)     8(r0st,x) 


3r     3r  . 

st     st 


3r, 


st 


3X 


(8-1) 


We  now  apply  the  Runge-Kutta-Gill  integration  formulas  to 
determine  9 (r  ,)/3(a)  from  its  time  derivatives,  given  by 


(7-5)  and  (7-11) 
K   = 


t  &   Y 
dt   p 


R  x   =  $  (K   -  6   Q  ) 


P    P 


> 


Y  _,n  =  Y  +  R  .- 

p+1  p    p+1 

Q  4.n  =  Q  +  3R  ,.  +  e   K 

XD  +  1  XD        D+l       D 


•p  +  1 


P  +  1       P    P  J 


P  =  0,1,2,3 


(8-2) 


After  each  cycle  an  intermediate  value  for  d/dt[3Cr  , )/3(a)] 

J  st 

is  computed  from  (7-5)  and  (7-11)  using  Y   for  3(r  .  )/3(.a) 

°   p        st 

(i.e.  for  3(r  ,  )/3(r0  .)  in  (7-5)  and  3  (r  J.)/3(x>  in  (7-11)). 
st      ust  st     A 

The  value  of  3  (r  . )/3(a)  for  t  +  At  is  then  Yu .   The  initial 

st  H 

value  for  Q   =  0o  is  a  6x (6+104)  zero-matrix  and  for  Y   =  Y0 

p     u  p     ° 

the  analytical  derivatives  at  the  starting  time  t0  as  given 
in  section  7. 
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5  , , , 5  , e   are 
p+1'  p'  p 


p 

bp 

6 

P 

EP 

0 

1/2 

2 

-1/2 

1 

1    -    A/2 

1 

A/2    J    1 

2 

1    +    A/2 

1 

-1    -    A/2 

3 

1/6 

2 

-1/2 

(8-3) 


This  method  is  also  used  for  the  computation  of  the 
starting  values  for  the  components  of  the  velocity  and 
position-vectors  with  Y0  now  as  the  initial  state  vector  as 
obtained  for  the  different  satellites  from  the  Smithsonian 
Institution.   Q0  again  is  a  zero  vector.   For  both  computations 
a  time  step  At  of  0.1  minute  is  used  only  in  the  starting  procedure 


2   The  Adams-Bashforth  Method 
Here  we  have 

t 


x 


L*-k+l 


'P+1 

=  /     P(t)dt 

t 
P 


=  At 


X 

E      Y      V 

y 

m=0 

•• 

LzJ 

(8-4) 


where  the  constants  y   are  independent  and  will  be  calculated 

m         r 

numerically  below.   At  is  the  time  step;  q  is  the  number 
of  steps,  which  in  this  case  is  10;  and  P(t)  is  the 


interpolating  polynomial  over  the  interval  [t    ,t  ]. 
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From  the  Runge-Kutta-Gill  method  initial  values  for  the 
state  vector  are  supplied.   Using  these  values 


X 

X 

X 

y 

5 

y 

5 

y 

.z_ 

_z_ 

_z 

Jq-i 


q-2 


x 

y 
l_2J 


are  computed 


from  the  equation  of  motion  so  that  the  differences  V 


m 


x 

y 

LZ, 


can  be  calculated  from  the  relations 


,m 


X 

X 

X 

y 

= 

y 

- 

y 

_z_ 

p 

_z_ 

P 

_z_ 

X 

/ 

X 

y 

= 

v| 

V-1 

y 

z 

p 

\ 

z 

p-1 


(3-5) 


The  expression  on  the  right  hand  side  of  equation  (S-M-)  is 


thus  known,  and 


can  be  calculated  for  p  >  q.   After 


p+1 


using  the  corrector  formula,  explained  in  sec.  8.3,  the  index  p 

is  increased  by  1,  and  the  same  formula  is  used  to  calculate 

x 

y      ,  etc.   With  the  help  of  the  following  recurrence 

L4J  p+2 

formula,  which  is  derived  in  Henrici  [1962,  p.  193],  numerical 


values  for  y   were  calculated 
'm 

i      i  : 

Ym  +  2  Ym-1  +  3  Ym-2  +  *  *  *  +  m+1  'm-q 


YT 


=  1,  m  =  0,1,2,  . 
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8.3   The  Adams -Moulton  Method 

Whereas  the  Adams -Bashforth  method  is  a  predictor  formula 
the  Adams -Moulton  method  is  a  corrector  formula.   Here  we 


have 

x 

y 


,--jp+i      "- 


p  +  l  q 

=  J  +     P(t)dt  =  At  I      yZ    V 


t 


m=0 


m 


"P 


x 

y 

L_Z_I 


(8-6) 


p  +  l 


where  P(t)  is  the  interpolating  polynomial  over  the  interval 

A 

[t    ,.  .  t  ,,].   The  recurrence  formula  for  y   is  similar  to 
p-q  +  1    p+l  'm 

that  one  for  y      and  can  also  be  found  in  Henrici  [1962]  on 
m 

page  194.   Formula  (8-6)  is  used  like  the  Adams-Bashforth 


formula  (8-4)  except  that  now  only  the  values 

are  known.   It  is  employed  to  redetermine 


x 

y 

z 

L.  _jp 


p-q  +  1 
an   iterative  procedure,  where  an  approximate  value 


x 

y 

* 

LJp  +  1 


p-1 
in 


x 


has  been  obtained  using  the  results  of  the  predictor  formulas 

(Adams-Bashforth  and  Cowell-Stormer ) 

(0) 
x 

-,m 


(1) 


x 

y 

L*Jp+l 


Calculating 


P 
(1) 

p+l 


+  At   E   y   v 

m=0 


X 

y 

l-Z_ 


(8-7) 


p  +  l 


and  re-evaluating  the  differences,  a 


better  values  is  then  obtained  for 


p+l 
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8 . 4   The  Stormer  Method 

By  integrating  a  differential  equation  of  the  second 
order  y"(x)  =  f(x,y,(x))  twice,  we  obtain  finally  the  formulas 
for  the  computation  of  the  position  vector 


-  2 


p+1 


L  Jp-1 


=  At2   E  a      V 

n   m 

m=  0 


x 

y 

L-Z_ 


(8-8) 


where  again  the  above  used  symbols  are  defined  as  before. 

a   is  obtained  with  the  help  of  the  following  recurrence 
m  to 


formula 

2  2 

0      =  1  -  x  h„  a  ..    -   Tr-  h0a   0  - 
m        3   2   m-1    4   3  m-2 

with  h   =l+i+...+i;m=l,2 

m        2  m         ' 


2 


h  . , an ;m=l , 2 


m+2   m+1  °' 


Formula  C8-8)  is  used  in  the  same  manner  as  the  Adams-Bashforth 
method. 

8.5   The  Cowell  Method 
Here  we  have 


X 

X 

X 

y 

-    2 

y 

+ 

y 

z 

z 

-r-,          -\ 

z 

p-1 


=  At 


p-2 


1  a"    Vm 
m=0   m 


where  a 


m 


2  .'•      2     * 

3  2   m-1    4   3  m-2 


X 

y 

2 

P 

2 
m+2 

-  h 
m 

(8-9) 


+  1  cr0 ,  m  =  1,2, 
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8.6   Integrating  Procedure  for  Variational  Equations 

The  system  of  variational  equations  (7-5),  (7-11)  can 
be  written  as 

R  =  SR  +  T 


where  R  = 


8r  ,  3r  . 

st    st 


(8-10) 
(8-11) 


and  is  a  6x(6+104)  matrix; 
3r 


S  = 


st 


8? 


,  and  is  6x6 


st 


and  T 


dt 


0 


st 


3X 


and  is  also  6*(6+104). 

The  equations  (8-10)  are  integrated  using  the  qth  order 

Adams-Moulton  formula.   At  the  (p+1)    step, 

r         st  r 

q        *       • 

R    xn    =    R      +    At    E      8        R    ... 
P+l  P  p-o      qp     P+l-P 


(8-12) 


where      3         =    (-1) 

qp 


C>>    Y*    ♦    C%h    Y*p+1    + 


+  <3>  y! 


p    =    0,1    . . . ,    q. 

q    =    0,1    ...     . 

Assuming  that  R  and  R  up  to  the  pth  step  are  known  from 
previous  applications  of  this  procedure  (or  at  the  start 
from  initial  values  obtained  by  the  Runge-Kutta-Gill  method), 
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we  see  that  the  only  quantity  not  known  on  the  right-hand 

side  is  R  ,,.   Let  us  rewrite  (8-12)  as 
p  +  1 

q 

R  ..  =  R   +  At  E   3*    R  xn    +  At  3*   R  in 
p  +  1     p        1   qp    p  +  l-p        qo   P  +  1 

Substituting  for  R  ...  from  (8-10),  we  obtain 

°      p+1 

q 

R  .-  =  R   +  At  t      3*  R  xn    +  At  3*   SR  xn+At  3*   T   , 
p  +  1     p  qp  p  +  l-p        qo    p  +  1      qo 

and  solving  for  R  ,,  , 
&      p  +  1 

R   xn    =    [i-At    3"    s]        (  R   +At  ("3"      T   +      £      3*      R   x.      1) 
p  +  1         L  qo   J        I     p  L     qo  p=1      qp      p  +  l-pj[ 


The  matrix  in  brackets  to  be  inverted  is  6x6,  and  the  matrix 

in  braces  is  6x(6+10U).   Following  the  computation  of  R  .-,,  ^D  +  i 

can  be  obtained  directly  from  (8-10)  for  use  in  the  next 

step. 

The  actual  computational  procedure  breaks  the  matrix  R 
into  components  as  shown  in  (8-11)  and  follows  the  outline 
given  by  Riley  et  al .  [1967]. 

8.7   The  lagrangian  Interpolation  Method. 

To  interpolate  between  the  time  steps,  Lagrange's 
interpolation  for  equidistant  abscissas  as  given  by  Henrici 
[1964]  on  pp.  201  -  202  is  applied.   Because  an  equal  time 
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step  At  is  used,  we  can  assume  that  the  points,  where  the 
values  (f,  )  of  the  positions,  velocities  and  variational 
equations  are  given,  are  equally  spaced.   Here  we  have 

n 
p(s)  =   E   £  (s)  f 

k=m 

where 

n 

l    f  O  =   IT    s  9, 

VS ;     J1    k-q 

q=m     ^ 

q*k 

This  representation  of  the  Lagrangian  polynomial  p(s)  is 
independent  of  the  size  of  At.   The  functions  )L  (s),  which 
Henrici  calls  the  normalized  Lagrangian  interpolation  co- 
efficients, depend  only  on  s  (the  relative  location  of  the 
observation  time  t  with  respect  to  next  lowest  time  t,  )  and 
on  the  integers  m  and  n,  which  are  the  bounds  of  the  set  of 
interpolating  points.   It  turned  out  that  the  results  for 

n=4-  and  m  =  -3  were  sufficiently  accurate. 

(tj) 
-3    -2    -1    0    1    2    3    1 

Using  this  method  only  the  values  f,  for  eight  time  steps 
and  the  value  t,  need  be  saved.   The  value  t,  is  incremented 
by  At  for  each  step  in  the  integration  process,  and  the  time- 
sequenced  observations  are  processed  as  soon  as  the  time  of 
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the  next  observation  (t  ,  )  is  exceeded  or  equaled  by  the 

obs 

current  t, .   We  then  have 

s  =  1  -  (t,  -  t,  )At 
1     obs 

The  f,  are  stored  in  a  circular  fashion  so  that  the  latest 
values  replace  the  former  values  for  f_„,  and  we  need  only 
monitor  the  index  of  the  current  location  of  f_„  in  the  eight 
place  array. 


^ •   Forces  Besides  the  Earth's  Gravity  Field  Acting  Upon  a 
Satellite 

In  addition  to  the  perturbations  caused  by  the  varia- 
tions of  the  earth's  gravitational  field  as  observed  in  the 
previous  chapters,  satellite  orbits  will  also  be  perturbed 
appreciably  by  the  gravitational  attractions  of  the  sun  and 
moon,  the  radiation  pressure  of  the  sun,  and  the  drag  due 
to  the  atmosphere. 

9.1  Attraction  of  the  Sun  and  the  Moon 

All  the  effects  of  lunar  and  solar  gravitation  upon  an 
artificial  satellite  may  be  expressed  by  adding  additional 
terms  for  the  potential  field  through  which  the  satellite 
travels.   If  r~  denotes  the  position  of  the  sun  or  moon  and 
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M„  its  mass,  the  disturbing  function  may  be  written  [Brouwer 
and  Clemence,  1961,  p.  465] 

AF  =  GM2(    1      -    -   T-^-\  (9-1) 

\  |r2  -  r|     r2   / 

The  first  term  gives  the  potential  due  to  the  disturbing 

body.   The  second  term  arises  from  choosing  the  center  of 

mass  of  the  earth  rather  than  the  center  of  mass  of  the 

entire  system  of  bodies  as  the  origin  [Danby   1962]. 

If  we  differentiate  (9-1)  with  respect  to  X  ,  where  X   stands 

r        m'        m 

for  the  coordinates  x,y,z,  we  get 

||£  =  GM,/-   -  Xm  -    +  C  (  ;   1  ; i-  A      (9-2) 

3Xm      2\    |?2  -  ?|3     m  |?2  -  ?|3     r23   / 


w 


th 


ith  Tr,    -   r„  (£,£;,£;),  where  L5Co  lie  in  the  equator  wi 

C-,  pointing  toward  the  vernal  equinox.   If  — 5-  is  factored 

r 
2 
out   we    get 

GM  r 

U£  =    -~  [-5      +    C5      -    x    )      -      2    -    y]  (9-3) 

8X  3  m  m  m       1  ->  -*  1  3 

m       r2  |r«    -  r  I 


with 


3  „-v   ->  3 

P2  r       2  2r*r2      r      ~2 

=    [1    +    (-    )Z    -    -    (£-)]  ^       •  (9-4) 


■*  -*-i  3  r>„  r   r0      r-, 

rrt    -   r  2  2         2 
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->  ->■ 

If  (— )   - (  — )  <  10"  ,  equation  (9-4)  will  then  be 

r2      r  r2   r2   —        ^ 

developed  in  a  series  so  that  we  get 

P2  *    n         3    r,rN2         2r'r,2    ,    rNn     ,    15r,rN2         2r'r,2    ,rN-,2 


1  .  4  [(JL)^  -  -^  (=£)]  +  J£[(JL)' £  (JL)] 


Ip   _  pi3        2    r2      r  r2   r2      8   r2      r  r2   r2 

.  35  [(JL)2  .  ^2   jj.  ]3  +  315  t±2    _   £^2   _r ,-,4       (9.5) 

16    r2      r  r    r2       128    r2      r  r2   r2 

The  quantity  r-rJlr  r„)  is  the  product  of  the  normal  vectors 

pointing  to  the  satellite  and  perturbing  body.   The  scale  of 

the  perturbing  body's  orbit  is  given  by  a2  and  is  needed  in  a 

few  of  the  small  terms,  where  r2  =  a~  (1  -  e  cosE).   The  series 

(9-5)  is  then  applied  to  (9-3)  to  obtain  the  difference 

3 

r2 
-5+5. 


'm    m  i  -*■    ■*■  i  3 

-p  _  -p 
|r2    r  | 

3 
The  coefficient  GM/r2  can  be  deduced  from  Kepler's  third 

law.   Applying  this  law  to  the  orbit  of  the  moon  and  the  sun, 

where  the  subscript  E  indicates  the  earth,  M  the  moon,  and  S 

the  sun  we  find 

G(ME  +  MM>  =  n^ 

G(ME  +  MM)  =  n^(l  -  eMcosEM)"3  (9,6) 

with  n„,  eM,  EM  being  elements  of  the  moon's  orbit  where 

n  =  mean  motion,  e  =  eccentricity,  E  =  eccentric  anomaly.   r„ 

stands  for  the  radius  of  the  moon's  orbit.   Hence, 

GM  M 

V1  =    M— +V  "M    (1    -    eMcosEM)-3  (9-7) 

rM  E  M 
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For  the  sun 


GM 


MS       2  -3 

"T  =  MQ  +  Mr   nS  (1  "  eS  cosES)  ^  ; 
rs     b     h 


(9-8) 


n^,  e„,  Eqj  Po  are  the  corresponding  quantities  for  the  sun's 
orbit  referred  to  the  earth;  the  ratio  M^/ (M„  +  M^)  is  very 
nearly  1. 

In  order  to  get  the  inertial  coordinates  £_,£_,£,- 
(eq   (9-3))  for  the  moon  and  sun,  respectively,  3  rotations 
have  to  be  applied 


5- 


with 


P1C-e)P3(-n)P1C-i)P3(-u) 


(cosE  -  e)/(l-e  cosE) 


/l-e2  sinE/(l-e  cosE) 
0 


P1(-e)  = 


P3(-fi)  = 


10     0 
0    cose  -sine 
0    sine   cose 

cosfi  -   sinfi   0 

sinfi     cosfi   0 

0        0     1 


•■ 

1 

0 

0 

P]_(-i)    = 

0 

cos    i 

-sin    i 

0 

sin    i 

cos    i 

COSU) 

-sinaj 

0 

P3(-oo)    = 

sinoo 

COSU) 

0 

0 

0 

1 

(9-9) 
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£  =  23°  26'  44'.'84 

fi,oj,  i  are  given  for  the  moon  by 

O-  =  12°  06'  46V05  -  0? 0529 5 386 52T 
u      =  196°  43'  52V316  +  0 .  164  358  002  5T 
cos  iM  =  0. 995970322 
sin  iM  =  0. 089683648 
with  T  =  MJD  =  3328  2.0 

For  the  sun,  i„  =  0. 

Hence,  P,  (-i)  =  I,  and  P~(-ft)  and  P„(-co)  may  be  combined  into 

one  transformation  employing  co  =  Q~    +  U}„. 

Q      +  cus  =  Sg  =  282°04'  45'.'92  +  0°  0000470684T 

For  the  solution  of  equations  (9-7),  (9-8)  values  for  e  and  E 

are  also  needed: 

eg  =  0.01677194 

e„  =  0. 054900489 
M 

With  the  help  of  m  =  E  -  e  sin  E  values  for  E  are  computed  by 
iteration.   The  mean  anomaly  m  is  found  by 

m   =  358°  0'  2V42  +  0° 98 56002647T 

mM  =  215°  31'  53V26  +  13° 0649924490T 
The  orbital  elements  given  here  for  the  sun  and  the  moon  have 
been  computed  for  1950.0  from  values  given  in  Explanatory 
Supplement  to  the  Astronomical  Ephemeris  and  the  American 
Ephemeris  and  Nautical  Almanac  [I960].   Only  secular  terms 
have  been  considered  here,  i.e.,  the  orbits  of  the  sun  and 
moon  have  been  approximated  by  rotating  ellipses. 
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9.2  Direct  Solar  Radiation  Pressure 

The  force  exerted  by  sunlight  pressure  is  proportional 
to  the  solar  flux  and  to  the  area  of  the  satellite  as  pro- 
jected on  a  plane  perpendicular  to  the  direction  of  the  flux. 
According  to  Shapiro  [1963]  the  acceleration  is 

K  A      r 

Q    =  Cjcr-)  (b    —  (9-10) 

xsp    M      c  r0 
r      sa        S 


where  rq  is  the  vector  from  the  earth  to  the  sun,  as  in  section 
9.1  ,  I  is  known  as  the  solar  flux  in  the  vicinity  of  the 


earth,  c  is  the  speed  of  light,  (A/M   )  is  the  area-to-mass  ratio 

sa 

of  the  satellite,  and  KR  is  a  function  lC(r,r«)  whose  value 
is  0  or  2  depending  on  the  location  of  the  satellite  in  its 
orbit.   In  the  space  cylinder  formed  by  the  earth's  shadow 
there  is  a  sharp  cutoff  of  the  force  for  any  point  within  the 
cylinder.   Here  the  value  K   =  0  is  used.   For  the  other 
part  of  the  orbit  the  value  2,  which  expresses  the  reflection 
characteristics  of  the  satellite,  is  used. 

Since  the  earth's  orbit  is  elliptical,  the  solar  flux 
will  vary  during  the  year: 

a   9 
I  =  I0  (— )  (9-11) 

E 
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where  a   is  the  semimajor  axis  of  the  earth's  orbit,  and  r_ 
e  L 

is  the  earth-sun  distance.   The  solar  constant,  I0 ,  which 

represents  the  mean  rate  at  which  enery  is  received  at  a  point 

2 
on  the  earth,  is  approximately  2  cal/cm   sec.  and  remains  quite 

constant.   Because  the  arcs  used  are  not  longer  than  7  days 
equation  (9-11)  was  not  applied. 

To  decide  whether  the  satellite  is  in  sunlight  or  in 
shadow  the  inner  product  of  the  position  vectors  of  the  sat- 
ellite and  the  sun  is  computed.   £  is  the  angle  between  these 
two  position  vectors.   The  satellite  is  in  shadow,  if  the 
following  conditions  are  fulfilled  (see  also  fig.  3) 


1)   cos  L    - 


J 


of  the  earth's 
omputation  of 
f f ects 
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The  components  of  eq   (9-10)  are  transformed  to  the  £-.9£9j£;_ 

coordinate  system  as  defined  in  section  9.1,  i.e.,  the  same 

rotations  are  used  as  for  the  computations  of  the  attraction 

3f 


of  the  sun.   The  partial 


KT 


as  necessary  for  eq   (2-5)  is 


numerically  integrated  in  the  same  manner  as  the  variational 
equations . 


9.3   Atmospheric  Drag 

Most  discussions  of  air  drag  in  the  current  literature 
on  artificial  satellites  begin  with  the  equation  for  the 
magnitude  of  the  drag  force  (e.g.,  Shapiro   1963  ) 

■  ,2 


with 


Qnd  =  "  2  CDP  CRA 


•) 


sa 


D 


(~ — )  =  area-to-mass  ratio  of  the  satellite 
sa 

C~  =  dimensionless  drag  coefficient 

here  approximately  2  and  a  variable  in  the 
adjustment 

p  =  air  density 

|rn|  =  relative  velocity  of  the  satellite  with 


D 


respect  to  the  atmosphere 


"  x    +    coy 

r*Di 

•*• 
rD    = 

y   -   cox 

= 

yD 

z           -1 

LzDJ 

where  co  is  the  average  angular  velocity  of  the  earth 
and  |?D|  =  /(*g  +  y£  I    z^ 


(9-17) 


C9-18) 
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The  drag  resistance  is  then 


ind 


2  CD  p  (R— }  'rD 
sa 


(9-19) 


The  air  density  p  is  given  as  a  function  decreasing'  with 
height 

p  =  p0.exp  [-(h  -  8-105)Ap]  , 

h  is  the  height  of  the  satellite  in  meters  obtained  from 
the  following  formula 


(9-20) 


h  =   r 


1  - 


R   . 
ma  3 


Ar2  +  E   z2) 
s 


(9-21) 


with  F   .  =    semimaior  axis 
ma  j  J 

2       2 
R   .  -  R  . 

and  E   =  _I3^L SiH  , 

R2. 
mm 


R  .   =  minor  axis  , 
mm  ' 

p0  in  equation  (9-20)  is  a  value  for  the  air  density  at  a 

height  of  80  0km  and  is  taken  from  the  U.S.  Standard  Atmosphere 

Supplements,  1966. 

5 
The  value  8-10   is  incorporated   because  of  the  chosen  p0 

and  the  unit  of  meters  for  h.   A   is  the  inverse  of  a  quantity 
known  as  the  "scale  height"  and  is  computed  by 


Ap  =  In  (-1)/(2.105) 
Po 


(9-22) 
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wh 


ere  p,  is  a  value  for  the  air  density  at  a  height  of  1000km 


and  is  also  taken  from  the  U.  S.  Standard  Atmosphere  Supple- 

5 
ments.   The  denominator  2*10   gives  the  height  difference 

between  p0  and  p,  in  meters. 

In  order  to  get  the  values  for  p0  and  p,  the  exospheric 

temperature  for  the  observation  time  has  to  be  determined. 

For  this  purpose  the  solar  flux,  the  monthly  solar  flux, 

and  the  geomagnetic  index  A   for  this  period  are  taken  from 

the  World  Data  Center  at  Boulder,  Colorado.   The  formulas 

used  here  can  be  found  in  the  U.  S.  Standard  Atmosphere 

Supplements.   First  the  variation  in  the  solar  cycle  is  taken 

into  account 

T0  =  362  +  3.60  F1Q  ?  (9-23) 

—  .  -22 

where  F, n  7  is  the  10.7cm  solar  flux  in  units  of  10 

2 
watts/m  /cycle/sec  averaged  over  three  solar  rotations.   The 

variation  within  one  solar  rotation  yields 

T'  =  T0  +  1.8  F10_?  -  F10_7)  C9-24) 

where  F, n  7  is  the  daily  solar  flux  averaged  over  the  ob- 
servation time  of  one  week.   Now  the  semiannual  variation  is 
supplied 

T0  =  TJ  +  f(d)  F  Q  (9-25) 
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where 

f(d)  =  (0.37  +  0.14  sin(2^  d  36551))  sin  Utt^^-) 

with  d  =  number  of  days  elapsed  since  January  1  of  each  year. 

The  correction  for  the  diurnal  variation  was  also 
supplied 

T  =  1.1  x  T0 

This  expression  is  a  very  simplified  version  of  the  equation 
for  the  diurnal  variation  found  in  the  U.  S.  Standard  Atmosphere 
Supplements . 

Finally  the  variation  with  geomagnetic  activity  yields 

T=T+AT  (9-26) 

00 

where  AT  =  a   +  ICQ  Cl  -  exp  C-0.08  A  3) 
P  P 

The  geomagnetic  index  A   is  averaged  over  the  observation 
time. 

With  the  value  for  the  exospheric  temperature  T  determined 
from  equation  (9-26),  values  for  the  air  densities  p0  and  p., 
are  taken  out  of  the  tables  6.1,  6.2,  6.3  given  in  the  U.  S. 
Standard  Atmosphere  Supplements. 

Values  for  the  air  density  computed  with  equation  (9-20) 
are  compared  with  values  given  by  Jacchia  [19  70]  and  found 
to  be  acceptably  in  agreement  over  the  range  of  the  height 
of  a  satellite  orbit. 
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In  order  to  obtain  the  partial  tttt-  in  eq   (2-5)  the 

9CD 

partial  derivative  of  Q  ,  with  respect  to  C-  is  numerically 
integrated  in  the  same  manner  as  the  variational  equations. 


1 0  .   Adjustment  of  Observations 

Rewriting  equation  (2-5)  and  using  a  vector-matrix 

notation  we  get 

A?  =  A  Ae  +  B  Ay  +  D  Ar   +  E  Ao   +  F  Ab  (10-1) 

^       s 


with 


Ae  = 


Ax 

Ay 

Az 
Ax 
Ay 
Az 


*■    H 


k  =  1,6 


ik 


Ae  being  a  (6x1)  vector  and  A  a  matrix  of  Cg  x  6)  elements, 
where  g  indicates  the  total  number  of  observation  equations 
per  arc . 

^Axi 

I   =  1,104 


AX 


A^ 


AX 


1  o  «t  J 


B=   i 


i£ 


"The  program  provides  an  option  to  do  the  adjustment  with  the 
parameters  for  air  drag  and  radiation  pressure  CCp.  and  KR ) 
as  constants.   In  that  case  equation  (10-1)  is  solved 
without  Ao. 
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with  B  being  a  matrix  of  (g  x  104)  elements 


Ar  = 
s 


Ar 


si  1 


Ar 


Sl  2 


Ar 


Sl  3 


Ar 


s2  1 


Ar 


Ar 


Ar 


ni 


n2 


n3 


Ax 


si 


Aysl 


Az 


s  1 


Ax 


S2 


Ax 


Ay< 


Az 


n 


n 


n 


D  = 


9f 


ia 


a  =  l,3n 


D  being  a  matrix  of  (g  x  3n)  elements  with  a=3(q-l)+m 
where,  as  explained  on  page  5,  n  stands  for  the  total  number 
of  stations,  m  denotes  the  mth  coordinate  of  r  ,  and  q 
indicates  the  qth  station.   Note  that  this  matrix  has  blocks 
of  zeros. 


Ao  = 


AK 


R 


AC 


D 


E  = 


V3KR)   >   \3CD) 


with  E  a  (g  x  2)  matrix 


Ab  = 


4bx 

• 

• 

Ab. 

•    1 

Ab, 

*■  n 


j  =  l,t 


ID 


Ab  being  a  vector  of  t  elements,  where  t  stands  for  the  total 
number  of  passes  per  arc  and  Fa  (g  x  t)  matrix.   Note  that  F 
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has  blocks  of  zeros. 

Equation  (lQ-1)  may  be  written  as 


[F,  G,  H] 


Ab 
Av 


=  t 


+   v 


C10-2) 


Here  I  is  a  vector  equivalent  to  Af  in  equation  (10-1)  and 
v  the  vector  of  the  residuals  where 


G  =  [A,E]  ,   H  =  [D,B] 


Av  = 


Ao 


,  and  Ari  = 


Ar 


By  means  of  the  covariance  matrix  S.  associated  with  the 


I 


observations  we  get  the  normal  equations 

Ab 


T   -1      T   -1      T   -1 
F  Z£  XF    F  Z£   G    F  E£   H 


T   -1      T   -1      T   -1 
G  Z£   F    G  E£   G    G  E£  ■LH 


T   -1      T   -1      T   -1 
HE.   F    H  E£  XG         H  E£   H 


Av 


An 


T   -1 

f  e  a 


T   -1 

G  h  * 


T   -1 
H  E£  I 


(10-3) 


In  the  formation  and  solution  of  the  least  squares 
normal  equations,  the  unknowns  fall  into  three  groups.   The 
three  groups  are:   bias  parameters  (one  Ab  for  each  pass), 
orbital  parameters  (orbital  elements  e,  plus  air  drag  C^, 
and  radiation  pressure  iC,),  and  surface  parameters  (station 

K 
coordinates  r   and  density  values  x)« 
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The  normal  equations  formed  from  all  observation 
aquations  having  the  form  of  (2-5)  contain  all  three  groups  of 
unknowns;  a  bias  parameter  for  every  pass,  a  set  of  orbital 
parameters  for  every  orbit,  and  one  set  of  station  coordinates 
and  densities.   Rather  than  directly  form  and  solve  these 
large  normal  equations,  we  eliminate  the  bias  and  orbital 
parameters  at  the  earliest  possible  stage.   Only  the  reduced 
normal  equations  associated  with  the  station  coordinates 
and  densities,  the  parameters  of  main  interest,  are  accumulated 
over  all  orbits. 

The  elimination  of  the  bias  and  orbital  parameters 
is  accomplished  by  an  application  of  Gaussian  elimination 
that  takes  account  of  the  special  structure  of  the  normal 
equations.   Generally,  the  elimination  of  X,  fro: 


>m 


M 


T 
P     Q 


~xl~ 

3   . 

u. 


u 


2  J 


(10-4) 
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produces  reduced  normals  for  X9  of  the  form 

(Q  -  PTM_1P)  X2  =  U2  -  PTM~1U1  . 


(10-5) 


If  desired,  X,  can  be  recovered  from  a  "back  solution" 


X]_  =  M  1(U1  -  P  X2)  . 


(10-6) 


If  M  is  "block  diagonal",  that  is,  of  the  form 


M  = 


m 
0 
0 


11 


0 


m 


22 
0     m 


33 


m 


nn 


(10-7) 


with  a  corresponding  partitioning  of  P  and  U,  , 


P  = 


m 


m 


m 


lc 
2c 
3c 


m 


nc 


U. 


u- 


u, 


u 


n 


the  reduced  normals  become; 
n 


Q  -   E   m.   m. .  m. 
j  =  1   DC   ]]   ]c 


n 


X, 


T    -1 
U0  -  l      m. .  m. .  u. 

2   .=1  11      ::  ] 


(10-8) 


Clearly,  each  term  of  the  summation  can  be  accumulated  as 
soon  as  m. • ,  m.  ,and  u.  are  available;  the  terms  can  be 


accumulated  in  any  order. 
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T   -1 
In  equation  (10-3),  F   Z0   F,  the  submatrix  of  the  normal 


■*•        T  -1 
equation  associated  with  the  bias  parameters  Ab,  and  G  Z0  G, 


•*■  T 

I 

the  submatrix  of  the  normal  equations  associated  with  the 
orbital  parameters  Av,  has  a  block  diagonal  form  which  is 
retained  after  the  elimination  of  Ab.   The  scheme  described 
above  is  applied  in  succession  to  eliminate  first  the  bias 
parameters  and  then  the  orbital  parameters.   The  contribution 
to  the  reduced  normals  arising  from  the  elimination  of  the  jth 
bias  parameter  is  computed  during  the  formation  of  the  partial 
normals  arising  from  the  jth  pass. 

After  all  observations  involving  the  1st  orbit  have  been 
processed,  the  accumulated  reduced  normal  equations  (reduced  by 
the  elimination  of  all  bias  parameters)  are  solved  for  pre- 
liminary values  of  corrections  to  the  orbital  parameters  of  the 
1st  orbit  (orbital  elements,  air  drag  and  radiation  pressure) 
and  surface  parameters  (station  coordinates  and  densities). 

The  normal  equations  are  then  reduced  a  second  time  to 
eliminate  the  orbital  parameters.   These  reduced  normal  equations, 
along  with  the  solution  for  the  densities  and  station  coordinates, 
are  stored  on  tape  for  introduction  during  the  formation  of 
normals  from  the  second  orbit  according  to  formula  (10-8). 

During  the  processing  of  the  observations  from  the  oth  orbit, 
the  reduced  partial  normals  from  this  orbit  (bias  parameters 
eliminated)  are  accumulated,  added  to  the  fully  reduced  normals 
carried  forward  from  the  (o-l)th  orbit,  and  solved  for  the 
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oth  preliminary  values  of  orbital  (for  the  oth  orbit)  and 
surface  parameters.   Then  the  oth  orbital  parameters  are 
eliminated  to  produce  reduced  normals  to  be  carried  forward, 
along  with  the  oth  preliminary  values,  to  the  (0+l)st 
solution.   The  final  solution  for  station  coordinates  and 
density  values  is  obtained  after  all  orbits  are  processed. 

In  fact,  the  accumulation  involved  in  the  formation  of 
reduced  normals,  as  above,  is  already  an  intrinsic  part  of 
most  Gaussian  elimination  programs,  such  as  the  Gauss-Jordan 
used  here.   Retrieval  of  the  reduced  normals  from  such  a 
program  requires  only  minor  modifications. 

The  results  of  the  adjustment  for  the  density  values  x5 

are  now  taken  to  compute  normalized  coefficients  C   and  S 

r  nm      nm 

in  order  to  compare  these  values  with  existing  satellite 

results.   If  C     and  S     denote  the  harmonic  coefficients 
nmu       nmu 

introduced  in  eq   (5-2)  to  define  the  normal  potential  U, 

we  obtain  [Koch   1968] 

n        104 

C    =  C     +  - £  y    \  \        n  — 

nm    nmu    fn    , .  s,  „   n  n    ,  AZ      {*      r   P   (sine}))  cos  mA  dE   (10-9) 
(2n+l)kM  a   £=1     AEp      nm 

-,  104 

S    =  S     + I    Xn      I!      rn  P   (sine}))  sin  mX   dE   (10-10) 

nmu    (2n+1)kM  an  ^H      ^ 

The  integral  over  the  surface  elements  AR  in  ((10-9)  and  (10-10)) 
is  solved  numerically  by  dividing  AE,  into  9  subdivisions. 
The  origin  of  the  earth-fixed  coordinate  system  and  of 
the  orbital  system--as  mentioned  earlier--is  the  center  of 
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mass  of  the  earth.   Hence  the  harmonic  coefficients  C,  n ,  C,  -.  , 
and  S  _ ,  must  equal  zero.   Furthermore  Co-i  and  S^-,  should  be 
small  in  comparison  to  the  rest  of  the  harmonic  coefficients 
since  during  the  orbit  computations  the  z-axis  effectively 
coincides  with  the  rotational  axis  of  the  earth.   To  insure 
this,  constraints  in  the  form  of  observation  equations  with 
small  variances  (see   Koch  and  Pope   1969)  are  set  up  ac- 
cording to  (eq  (10-9)  and  (10-10))  and  their  contribution  is 
added  to  the  normal  equations.   Five  such  constraint  equations 
are  used,  setting  each  of  C,  n ,  C,  ,-j  S-,  ,  ,  C„,  ,  and  S~-.  equal 

to  zero.   As  an  option,  C    can  also  be  constrained. 

r     '   oo 

In  addition  to  these  constraints  the  longitude  of  one 
station  was  held  fixed  in  the  solution  to  prevent  a  singularity 
which  would  arise  if  all  the  station  longitudes  and  the  right 
ascension  of  each  orbit  were  all  unknown.   We  get  the  follow- 
ing condition  equation 

j-        y    +    Ay  .         y 

arc    tan   - — t — tt~  -    arc    tan   2- 

x  +  Ax  x 


or  by  linearization 


&■  -  2$&  )   =  o  .  (lo-ii) 


1  ♦  C^)2  V   x   x2   / 

X      v 

For  some  of  the  stations,  which  have  been  moved  during  the 
observation  time  of  different  orbits,  the  surveyed  distances 
between  the  old  and  new  location  are  introduced  as  constraints. 
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The  differences  in  the  coordinates  (Au,Av,Aw)  are  weighted 
based  on  a  mean  square  root  error  according  to  the  obtained 
accuracy  and  observation  equations  with  these  differences 
are  added. 
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